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ABSTRACT 

A  quantum-mechanically-based  potential  energy 
function  that  describes  interactions  of  dimers  of  the 
explosive  cyclotrimethylenetrinitramine  (RDX)  is  used  to 
predict  polymorphic  structures  of  crystalline  RDX.  The 
potential  was  first  used  in  crystal  structure  prediction 
methods  which  generate  polymorphs  of  RDX  and  provide 
a  0  K  ordering  in  energy.  The  13  low  energy-structures 
generated  by  crystal  structure  prediction  methods  were 
subsequently  subjected  to  isothermal-isostress  molecular 
dynamics  (N.sT  MD)  simulations  to  establish  stability 
rankings  and  assess  effect  of  inclusion  of  temperature. 
For  both  crystal  structure  prediction  methods  and 
molecular  dynamics  simulations,  the  low-energy 
polymorph  corresponds  to  the  experimental  structure  with 
crystallographic  parameters  in  outstanding  agreement 
with  experiment. 

1.  INTRODUCTION 

Advances  in  quantum  mechanical  methods  and 
computational  resources  have  allowed  ARL  researchers  to 
develop  powerful  predictive  capabilities  for  use  in  the 
design  and  development  of  advanced  energetic  materials. 
These  computational  tools  are  used  to  predict  properties 
of  energetic  materials  related  to  performance  or  hazard, 
thus  allowing  for  assessment  of  a  notional  material  using 
computational  methods  only.  One  key  property  is  the 
crystal  density;  we  have  expended  considerable  effort  in 
developing  a  robust  capability  for  its  prediction  using 
theoretical  information  only. 

Although  the  molecular  structure  of  single  molecules 
can  be  easily  and  accurately  predicted  using  highly- 
accurate  quantum  mechanical  (QM)  methods,  the  same 
cannot  be  said  for  predicting  the  system  in  its  crystalline 
state.  In  fact,  prediction  of  a  crystal  structure  using 
information  only  about  the  atomic  arrangement  of  a 
constituent  molecule  was  once  considered  impossible 
[Dunitz,  2003]  due  to  substantial  theoretical  and 
computational  challenges,  including  properly  representing 


the  dominant  binding  forces  within  molecular  crystals, 
van  der  Waals  interactions.  While  these  forces  can  be 
modeled  ab  initio  by  wave-based  methods  of  quantum 
chemistry,  the  high  computational  cost  of  such 
approaches  has  prevented  extensive  applications  of  these 
methods  to  molecular  crystals.  The  only  computationally- 
feasible  quantum-mechanically-based  approach  to  treat 
molecular  crystals  is  Density  Functional  Theory  (DFT), 
shown  to  be  inaccurate  for  systems  for  which  van  der 
Waals  forces  are  dominant  (e.g.  molecular  crystals)  [Byrd 
et  al.,  2004]. 

These  computational  restrictions  have,  to  this  point, 
forced  us  to  include  empirical  models  in  methods  for 
prediction  of  crystal  structures.  These  empirical  models 
rely  on  experimental  information  for  parameterization  and 
are  limited  in  their  predictive  capability.  To  illustrate,  in 
a  previous  study,  we  used  a  molecular  simulation 
approach  that  allows  for  the  a  priori  prediction  of  the 
crystal  density  and  crystallographic  parameters  using 
information  about  a  single  molecule  (i.e.  ab  initio  crystal 
prediction);  however,  the  intermolecular  interactions  were 
described  by  this  empirical  potential  [Rice  and  Sorescu, 
2004].  A  total  of  174  CHNO  molecular  crystals  for 
which  experimental  crystallographic  information  is 
available  were  subjected  to  this  approach.  The 
calculations  produced  148  crystals  whose  crystallographic 
parameters  and  molecular  configurations  matched  those 
of  the  experimental  counterpart.  Unfortunately,  the 
method  and  potential  failed  to  identify  26  of  the 
structures.  Subsequent  application  of  the  method  and 
model  to  high-nitrogen  molecular  crystals  was  completely 
unsuccessful.  These  failures  demonstrate  that  atomistic 
simulations  using  empirical  models  often  produce 
unsatisfactory  results  in  describing  systems  outside  of 
those  included  in  the  fitting  data. 

Here  we  provide  an  alternative  approach  to  such 
empirical  modeling  and  demonstrate  it  for  the  well-known 
energetic  material,  cyclotrimethylenetrinitramine  (RDX). 
The  potential  energy  function  for  describing  crystalline 
RDX  is  based  on  dimer  interactions  calculated  using  a 
highly  accurate  quantum  mechanical  method,  symmetry- 
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adapted  perturbation  theory  based  on  the  density- 
functional  description  of  the  monomers  [SAPT(DFT)] 
[Podeswa  et  al,  2007].  In  order  to  establish  the  quality  of 
the  SAPT(DFT)-based  potential  and  its  utility  for  use  in 
molecular  simulations  of  condensed  phase  RDX,  we 
subjected  it  to  a  series  of  increasingly  rigorous  atomistic 
simulations  to  determine  its  suitability  for  use,  its 
deficiencies  and  its  limitations.  As  a  first  step  in  such  an 
assessment,  we  used  it  to  model  crystalline  RDX  at 
ambient  conditions  using  isothermal-isostress  molecular 
dynamics  (N.s'T-MD)  simulations  assuming  as  an  initial 
configuration  the  structure  of  the  experimental  crystal. 
Subsequent  tests  included  generation  and  energy  ranking 
of  polymorphs,  and  subsequent  subjection  of  the 
polymorphs  to  N.s'T-MD  simulations  for  energy  ranking 
when  thermal  motion  is  included. 


2.  COMPUTATIONAL  DETAILS 

A  complete,  six-dimensional  potential  energy  surface 
(PES)  describing  interactions  between  dimers  of  RDX 
was  generated  using  over  1000  single-point  calculations 
using  the  SAPT(DFT)  method  [Misquitta  et  al,  2005].  A 
functional  form  of  the  potential  energy  surface  assumed  a 
site-site  formulate  over  all  nuclei  of  a  pair  of  dimers,  and 
was  fitted  to  the  SAPT(DFT)  points.  Details  of  the 
calculations  and  the  fitting  are  given  in  [Podeswa  et  al, 
2007].  Before  subjecting  the  PES  to  polymorphic  crystal 
structure  prediction,  we  used  it  in  isothermal-isostress 
molecular  dynamics  (N.vT-MD)  simulations  using  the 
experimentally  determined  crystal  as  the  initial  state. 

Polymorphic  crystal  structures  corresponding  to  local 
minima  on  the  PES  were  identified  using  the  method  of 
ab  initio  crystal  prediction  as  implemented  by  Ammon 
and  co-workers  [Holden  et  al.,  1993]  in  the 
MOLPAK/WMIN  suite  of  crystal  prediction  software. 
The  MOLPAK/WMIN  procedure  is  similar  to  those  of 
most  ab  initio  crystal  prediction  methods,  which 
generally  consists  of  three  steps,  with  the  first  step 
corresponding  to  generating  a  three-dimensional  model  of 
the  molecule  that  will  be  used  to  construct  candidate 
crystals  of  different  symmetries  (denoted  hereafter  as  the 
“test  molecule”).  Although  quantum  mechanical 
procedures  can  reliably  and  rapidly  predict  geometries  of 
the  model  used  in  this  first  step,  the  three-dimensional 
model  we  used  corresponds  to  the  experimental  structure 
at  room  conditions  [Choi  and  Prince,  1972]  since  this  was 
the  configuration  used  for  the  generation  of  the  SAPT- 
DFT  points  to  which  the  PES  was  fitted. 

The  second  step  of  the  procedure  involves  creating 
hypothetical  crystal  structures  using  the  molecular 
models.  In  the  MOLPAK  procedure,  an  initial  packing 
arrangement  is  obtained  by  surrounding  the  test  molecule 
with  a  coordination  sphere  containing  other  molecules. 


The  contents  and  three-dimensional  structure  of  the 
sphere  are  dependent  on  the  crystalline  space  group 
symmetry.  The  definitions  of  the  various  coordination 
spheres  used  in  a  MOLPAK  calculation  were  obtained 
from  detailed  analyses  for  a  large  number  of  organic 
crystal  structures  [Holden  et  al.,  1993].  The  analyses 
showed  that  the  most  probable  number  of  molecules  in 
the  coordination  sphere  is  14,  and  that  specific  “patterns 
and  sub-patterns”  were  apparent  in  the  three-dimensional 
structure  of  the  coordination  spheres  (Holden  et  al., 
1993).  The  version  of  MOLPAK  used  in  this  study 
samples  51  coordination  groups  corresponding  to  the 
triclinic  (PI,  PI),  monoclinic  (P2b  P2!/c,  Cc,  C2,  C2/c, 
Pc,  P2/c,  P2i/m,  P2/m,  P2,  Pm,  P2/m),  and  orthorhombic 
(Z=4,  P2i2]2,  P2i2,2i,  Pca2i,  Pna2t,  Pnn2,  Pba2,  Pnc2, 
P222j,  Pmn2i,  Pma2,  Z=8,  Pbcn,  Pbca)  space  groups. 

Adjacent  molecules  in  the  initial  coordination  sphere 
are  then  systematically  moved  in  small  steps  toward  the 
centrally-located  test  molecule.  At  each  step,  the 
potential  energy  of  the  “crystal”  is  calculated.  This 
continues  until  a  repulsion  criterion  is  met.  Once  this 
criterion  is  met,  the  packing  procedure  stops,  and  the 
volume,  crystallographic  parameters  and  Eulerian  angles 
that  describe  the  orientation  of  the  test  molecule  are 
stored. 

Next,  a  new  candidate  crystal  is  generated  by  first 
reorienting  the  test  molecule  about  only  one  of  the 
Eulerian  axes  by  10°,  and  repeating  the  packing  procedure 
described  heretofore  until  the  entire  Eulerian  space  is 
sampled.  Rotations  in  10°  steps  about  the  three  Eulerian 
axes  will  result  in  the  generation  of  6,859  (193) 
orientations  and  hypothetical  crystal  structures  for  each  of 
the  possible  space  group/coordination  sphere 
combinations. 

After  the  -7000  structures  were  generated  for  each 
coordination  sphere  geometry,  they  are  ranked  according 
to  density.  The  500  most  dense  structures  are  subjected  to 
full  energy  minimization,  where  minimization  is 
performed  with  respect  to  the  crystallographic  parameters. 
The  energy  minimization  is  space-group  symmetry 
restricted  (i.e.,  the  space  group  symmetry  is  conserved 
throughout  the  minimization),  and  performed  using  the 
code  WMIN  [Busing,  1981]. 

At  the  completion  of  the  energy  minimizations  of  the 
500  most  dense  hypothetical  crystals  generated  for  each 
of  the  5 1  space  group/coordination-sphere  geometries,  the 
user  has  the  information  needed  to  identify  the  “correct” 
crystal  structure  according  to  his  specific  criteria  (e.g. 
lowest  lattice  energy  or  highest  density  [Lommerse  et  al., 
2000]).  For  our  purposes,  we  have  ranked  the  crystals 
according  to  lowest  lattice  energy. 

The  MOLPAK/WMIN  procedure  proceeds  under  the 
assumption  that  all  important  regions  of  configuration 


2 


space  are  adequately  sampled  and  that  the  interatomic 
interactions  are  sufficiently  accurate  to  describe  packing 
and  lattice  energies. 

Due  to  limitations  in  the  MOLPAK/WMIN  software, 
a  scaled  down  version  of  the  SAPT(DFT)  potential  was 
used  in  the  prediction  of  the  polymorphs;  its  form  is  a 
standard  pair  additive  exponential  six  form  with 
coulombic  interactions,  required  by  our  version  of  the 
MOLPAK/WMIN  software.  It  will  be  denoted  hereafter 
as  the  SAPT(DFT)exp_6  potential.  While  this  modified 
version  of  the  PES  was  utilized  in  the  polymorph 
predictions,  its  features  were  compared  with  those  of  the 
full  SAPT(DFT)  potential,  and  showed  only  minor 
differences  in  most  cases.  We  checked  this  by 
performing  one  isothermal-isostress  molecular  dynamics 
(N.S'T  MD)  simulation]  on  the  low-energy  polymorph  with 
both  potentials:  The  cell  vectors  of  the  time-averaged 
geometry  differed  by  no  more  than  0.7%  in  the  two  cases. 

The  13  lowest-energy  cells  generated  through  the 
WMIN  optimization  are  then  used  as  initial  structures  in 
N.yT-MD  simulations  to  determine  stability  and  energy 
ranking  relative  to  the  experimental  structure.  For  the 
N.yT-MD  simulations,  the  full  SAPT(DFT)  PES  was  used. 
The  rigid-molecule  approximation  was  assumed  and  no 
symmetry  constraints  were  imposed  during  the  trajectory 
integration.  The  DL  POLY  version  2  suite  of  molecular 
dynamics  software  was  used,  with  thermostat  and  barostat 
relaxation  times  set  at  2.0  and  5.0  ps,  respectively. 
Supercells  for  the  polymorphs  consisted  of  blocks  of  unit 
cells,  with  sizes  selected  to  ensure  that  the  perpendicular 
widths  between  opposing  faces  of  the  simulation  cell  were 
at  least  twice  the  interaction  potential  cutoff  distance 
(12.0  A).  Coulombic  interactions  were  handled  using 
Ewald  summations.  Atomic  velocities  corresponding  to 
T=300  K  were  then  assigned  and  a  trajectory  of  30,000 
steps  (1  time  step  =  0.001  fs)  was  performed  to  equilibrate 
the  system.  During  equilibration,  atomic  velocities  were 
scaled  at  every  fifth  step  to  more  quickly  drive  the  system 
to  the  target  temperature. 

After  equilibration,  time-averaged  thermodynamic 
properties  and  crystal  parameters  were  determined  from 
results  of  100,000-step  N.yT-MD  trajectory.  Also,  at 
every  500th  step  during  the  trajectory,  atomic 
configurations  were  recorded  to  generate  time-averaged 
information  about  the  contents  of  the  simulation  cell. 
Results  are  given  in  terms  of  center-of-mass  fractional 
(sx,sy,sz)  and  the  Euler  angles  (0.4>,v)  that  transform  the 
principal  axes  of  inertia  of  each  molecule  to  its  space- 
fixed  coordinate  system  (i.e.,  the  Cartesian  axes).  For 
each  polymorphic  form  studied,  the  molecular  parameters 
for  each  of  the  symmetry-equivalent  molecules  in  the  unit 
cell  were  averaged  over  time  and  over  all  unit  cells  within 
the  simulation  supercell. 


A  key  test  of  an  intermolecular  interaction  potential 
to  describe  molecular  crystals  it  is  ability  to  maintain  the 
proper  crystalline  space  group  symmetry  during  N.yT-MD 
simulations  at  the  appropriate  temperature  and  pressure. 
We  assessed  whether  such  was  maintained  during  the 
trajectories  by  generating  “ideal”  unit  cells  for  the  13 
polymorphs  subjected  to  N.yT-MD  simulation.  Each  ideal 
cell  was  constructed  using  the  time-averaged  atomic 
positions  of  one  of  the  molecules  in  the  simulation  cell, 
and  generating  the  appropriate  number  of  symmetry 
equivalents  required  to  construct  a  unit  cell  according  to 
the  space  group  symmetry  requirements.  The  molecular 
parameters  (center-of-mass  fractional  and  orientational 
Euler  angles)  of  the  time-averaged  unit  cells  generated 
from  the  N.yT-MD  simulations  were  then  compared  with 
those  of  the  ideal  crystals. 

3.  RESULTS  AND  DISCUSSION 
3.1  SAPT(DFT)-based  PES 

In  the  SAPT(DFT)  method,  the  interaction  energy 
can  be  decomposed  into  components,  each  of  which  has  a 
clear  physical  interpretation  and  describe  various 
interactions,  including  electrostatic,  exchange,  induction, 
exchange-induction,  dispersion,  and  exchange-dispersion 
interactions.  These  calculations  represent  the  largest 
system  ever  investigated  with  any  high-level  electronic 
structure  method.  5 1  local  minima  were  identified  on  the 
PES  for  this  system;  each  was  analyzed  in  terms  of  the 
physical  components  of  the  interaction  energy.  As 
expected,  the  dispersion  represented  the  dominant 
component,  although  in  some  cases,  the  electrostatic 
interaction  is  significant.  The  analysis  of  the  SAPT(DFT) 
PES  indicates  that  the  neglect  of  the  dispersion 
contribution  would  result  in  excessively  weakly-bound  (or 
unbound)  dimers.  This  analysis  provides  a  definitive 
explanation  as  to  why  conventional  DFT  treatments 
(which  do  not  properly  treat  dispersion)  provide 
inaccurate  descriptions  of  crystalline  RDX. 

As  a  first  step  in  assessing  the  utility  of  this  PES  for 
use  in  molecular  simulations  of  crystalline  RDX,  an  N,yT- 
MD  simulation  of  the  experimental  crystal  structure  at 
room  conditions  was  performed.  The  resulting 
crystallographic  parameters  and  orientations  of  the 
molecules  within  the  unit  cell  were  in  excellent  agreement 
with  experiment.  Figure  1  shows  the  superposition  of 
the  predicted  unit  cell  onto  the  experimental  one.  The 
agreement  with  experiment  is  excellent,  even  better  than 
typically  achieved  in  simulations  with  empirical  potentials 
fitted  to  a  given  set  of  crystal  data. 
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Figure  1.  The  predicted  RDX  unit  cell  superimposed  onto 
the  experimental  structure. 

3.2.  Prediction  of  RDX  Polymorphs 

MOLPAK/WMIN  predictions  of  the  density  and 
lattice  energy  (per  molecule)  for  the  low  energy  crystal 
for  each  of  the  51  coordination  groups  sampled  in  this 
study  were  performed,  with  the  lowest  energy  crystals 
corresponding  to  the  Pbca,  P2J/C,  P212l21,  Pna2h  and 
C2/c  crystalline  space  groups.  The  low  energy  crystal, 
corresponding  to  Pbca  (the  experimental  space  group),  is 
lower  in  energy  by  at  least  1.94  kcal/mol  than  all  other 
crystals,  thus  eliminating  any  ambiguity  in  the  polymorph 
ranking.  Of  the  thirteen  polymorphs  identified  by  the 
MOLPAK/WMIN  procedure  and  subsequently  subjected 
to  NsT-MD  simulation,  one  produced  an  amorphous 
structure  upon  introducing  thermal  and  dynamic  motion; 
this  corresponded  to  MOLPAK/WMIN  Solution  444, 
coordination  group  AK  (space  group  P2j/c).  This 
indicates  that  the  MOLPAK/WMIN  SAPT(DFT)exp_6 
structure  is  not  a  stable  minimum  on  the  SAPT(DFT) 
PES.  No  further  analysis  of  this  structure  was  performed. 

Deviations  of  the  NsT-MD  results  from  those  of  the 
MOLPAK/WMIN  predictions  were  no  greater  than  5.6% 
in  cell  dimension,  1.3°  in  cell  angle,  and  8%  in  density.  A 
more  rigorous  comparison  of  the  contents  of  the  cells 
generated  through  NjT-MD  simulations  and 
MOLPAK/WMIN  calculations  are  performed  through 
superposition  of  the  N.sT-MD  time  averaged  unit  cells 
onto  the  MOLPAK/WMIN  counterpart  as  described  in 
Rice  and  Sorescu  [2004].  The  maximum  rigid-body 
rotational  displacement  and  translational  displacement  of 
any  of  the  symmetry-equivalent  molecules  in  the  time- 
averaged  unit  cell  relative  to  those  of  the 
MOLPAK/WMIN  counterpart  were  16.3°  and  0.58  A, 


respectively.  However,  these  values  corresponded  to 
higher-energy  polymorphs.  Among  the  four  lowest- 
energy  polymorphs,  the  maximum  rigid-body  rotational 
displacement  and  translational  displacement  of  any  of  the 
symmetry-equivalent  molecules  in  the  time-averaged  unit 
cell  relative  to  those  of  the  MOLPAK/WMIN  counterpart 
are  8.7°  and  0.30  A,  respectively. 

The  N.vT-MD  lattice  energies  are  all  higher  than  the 
MOLPAK/WMIN  results  but  the  densities  are 
consistently  lower.  These  trends  can  be  partially 

attributed  to  thermal  effects  introduced  in  the  MD 

simulations,  whereas  the  MOLPAK/WMIN  calculations 
represent  a  0  K  result;  thermal  expansion  should  generate 
lower  densities  and  higher  lattice  energies.  Also,  there 
are  slight  differences  in  absolute  energies  between  the 
SAPT(DFT)exp-6  and  SAPT(DFT)  potentials,  resulting  in 
differing  relative  energy  rankings  of  the 

MOLPAK/WMIN  and  N.vT-MD  structures.  However, 

for  both  sets  of  calculation,  the  low-energy  structure 
corresponds  to  the  experimental  crystal. 

Table  1  provides  the  time-averaged  crystallographic 
parameters  for  the  4  lowest-energy  crystals  produced 
from  N.vT-MD  simulations  using  the  SAPT(DFT) 
potential. 


Table  1:  Crystallographic  parameters  calculated  using 
N  sT -MD/SAPT  (DFT)  methods.  _ _ 


Solution/Space 

Group 

a  (A) 

b  (A) 

c(A) 

Experiment 

13.82 

11.574 

10.709 

CB-4/  Pbca 

13.234 

11.609 

10.720 

CB17/Pbca 

13.377 

11.682 

10.755 

AM17/P2!/c 

8.132 

10.702 

9.662 

AM75/P2!/c 

7.236 

11.683 

10.719 

Solution/Space 

Group 

a(°) 

PO 

r(°) 

Experiment 

90.00 

90.00 

90.00 

CB-4/  Pbca 

90.00 

89.99 

90.00 

CB17/Pbca 

90.00 

89.99 

90.00 

AM17/P2i/c 

90.00 

96.41 

89.99 

AM75/P2!/c 

90.00 

71.46 

90.00 

The  initial  configurations  for  the  MD  simulations  were 
generated  from  calculations  using  the  MOLPAK/WMIN 
procedure  and  the  SAPT(DFT)exp_6  potential.  Table  2 
provides  the  maximum  deviations  of  the  monomer’s 
positions  and  orientations  in  the  time-averaged  unit  cells 
from  those  of  the  “ideal”  RDX  crystals  for  each  of  the 
four  low-energy  polymorphs.  Of  the  13  crystals  subjected 
to  N.vT-MD  simulations,  the  maximum  deviation  of 
center-of-mass  fractional  coordinates  from  ideality  among 
all  of  the  polymorphs  is  0.0032  fractional  units  (Solution 
75,  AM)  and  the  maximum  deviation  of  Euler  angles  from 
ideality  is  3.01°  (Solution  304,  AK).  The  good  agreement 
of  the  N.vT-MD  averaged  values  with  those  of  “ideal” 
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crystals  demonstrate  that  the  crystalline  space  group 
symmetries  are  well  maintained  throughout  the 
simulations  for  these  polymorphs.  Values  in  Tables  1 
and  2  are  given  in  order  of  increasing  lattice  energy,  with 
crystal  CB-4/Pbca  having  the  lowest  lattice  energy. 


Table  2:  Maximum  absolute  deviations  of 
monomers’  positions  and  orientations  from 
corresponding  values  in  an  ideal  RDX  crystal. _ 


Solution/Space 

Group 

Hv.vVIax 

AYjMax 

ASZ  Max 

CB-4/  Pbca 

3.1e-04 

2.5e-04 

1.9e-04 

CB17/Pbca 

2.3e-04 

5.2e-04 

3.0e-04 

AM17/P2i/c 

6.8e-04 

2.5e-04 

2.4e-04 

AM75/P2!/c 

3.2e-03 

1.8e-04 

7.8e-04 

Solution/Space 

Group 

A©  Max  (°) 

•dt^viax  (°) 

a'Fm,,  o 

CB-4/  Pbca 

0.16 

0.07 

0.13 

CB17/Pbca 

0.13 

0.14 

0.12 

AM17/P2i/c 

0.19 

0.64 

0.97 

AM75/P2!/c 

0.34 

0.02 

0.97 

Lattice  parameters  for  the  experimental  crystal  are 
also  given  in  Table  1  for  comparison  with  the  low-energy 
time-averaged  structure  (Solution  CB-4)  generated  in  the 

N. yT-MD  simulations.  Note  that  the  predicted  low- 
energy  structure  has  the  same  space  group  as  the 
experimental  crystal.  Agreement  with  experiment  for  all 
cell  parameters  is  very  good.  The  deviations  from 
experiment  of  the  N.yT-MD  densities  and  cell  edge  lengths 
in  the  a-,  b-  and  c-  directions  are  -0.8%,  0.4%,  0.3%,  and 

O. 1%,  respectively.  A  comparison  of  time-averaged 
center-of-mass  fractional  coordinates  and  orientational 
Euler  angles  for  the  eight  symmetry  equivalent  molecules 
in  the  unit  cell  with  corresponding  experimental  values 
for  RDX  at  298  K,  1  atm,  showed  outstanding  agreement. 
The  largest  deviation  in  molecular  orientation  occurs  for 
the  Euler  angle  vj/,  which  differs  from  experiment  by 
-2.1°.  Also,  the  largest  deviation  of  the  location  of  the 
molecular  mass  center  is  less  than  0.07  A. 

In  order  to  estimate  the  uncertainties  of  the 
SAPT(DFT)  PES,  we  evaluated  the  main  neglected 
component,  the  non-additive  three -body  contributions  to 
the  lattice  energy  for  the  experimental  geometry  of  the 
crystal,  which  is  the  main  component  of  the  energy  that 
was  neglected  in  these  calculations.  This  contains  the  E1F 
energy  (0.76  kcal/mol),  the  dispersion  energy  (1.00 
kcal/mol),  and  the  leading  asymptotic  dispersion  term 
(0.31  kcal/mol).  These  values  were  obtained  using  28,  7, 
and  ~  8  x  106  symmetry-unique  trimers.  Since  these  three- 
body  contributions  are  small  relative  to  the  two-body 
component,  the  lattice  parameters  should  change  very 
little  upon  their  inclusion  in  the  PES.  Additionally,  the 
total  three-body  effect  is  dominated  by  the  dispersion 
energy,  which  is  fairly  isotropic  and  would  be  similar  for 
all  of  the  polymorphs.  Thus,  we  conclude  that  the  relative 


energies  of  the  polymorphs  are  unlikely  to  change  after 
including  the  many-body  forces.  The  absolute  value  of  the 
lattice  energy  changes  by  6.7%  for  structure  CB/4, 
suggesting  that  the  good  agreement  of  the  pairwise 
additive  energy  with  experiment  is  due  to  a  cancellation 
of  errors,  including  basis  set  incompleteness  error  of  the 
pair  potential  countermanded  by  the  three-body  effects  of 
opposite  sign. 


CONCLUSIONS 

We  have  demonstrated  that  the  SAPT  (DFT)  method 
is  capable  of  producing  force  fields  for  interactions  of  the 
molecular  crystalline  explosive  RDX,  and  appears  to  be 
suitable  to  enable  reliable  predictions  of  crystal  structures 
of  similar  molecular  crystalline  compounds.  The 
predictions  were  done  entirely  from  first  principles,  with 
no  reliance  on  empiricism,  thus  demonstrating  that 
compounds  for  which  experimental  data  are  unavailable 
can  be  explored  using  theoretical  simulation 
methodologies.  We  expect  this  method  to  find  broad 
applications  not  just  in  energetic  materials  design  and 
development,  but  also  in  areas  of  crystal  design,  in 
particular,  to  screening  novel  materials  and  drug 
candidates,  screening  molecules  for  co-crystallization, 
and  identification  of  low-energy  polymorphs  of 
pharmaceutical  compounds. 
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